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Pulsar Emission above the Spectral Break - A Stacked Approach 
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NASA’s Fermi space telescope has provided us with a bountiful new population of gamma-ray sources following 
its discovery of over 150 new gamma-ray pulsars. One common feature exhibited by all of these pulsars is the 
form of their spectral energy distribution, which can be described by a power law followed by a spectral break 
occurring between and ~8 GeV. The common wisdom is that the break is followed by an exponential cutoff 
driven by radiation-reaction-limited curvature emission. The discovery of pulsed gamma rays from the Crab 
pulsar, the only pulsar so far detected at very high energies (E>100 GeV), contradicts this “cutoff” picture. 
Here we present a new stacked analysis with an average of 4.2 years of data on 115 pulsars published in the 
2nd Fermi-LAT catalog of pulsars. This analysis is sensitive to low-level ~100 GeV emission which cannot be 
resolved in individual pulsars but can be detected from an ensemble. 


1. Introduction 

One common feature exhibited by all known 
gamma-ray pulsars is the form of their spectral en¬ 
ergy distribution (SED) which can be described by a 
power-law followed by a spectral break occurring be¬ 
tween 1 and 8 GeV [Abdo et al.| 2013 . The unanim¬ 
ity of the break energy across the entire Fermi-hAT 
pulsar sample is suggestive that the sites of acceler¬ 
ation and processes of gamma-ray emission are com¬ 
mon across different pulsar types and that they are 
not strongly dependent on the pulsar spin or ener¬ 
getics. Further, it has been shown that across the 
Fermi-UAT pulsar sample the spectral-break energy 
is weakly correlated with the magnetic-field strength 
at the light cylinder Abdo et al. 1 2010b 2013 . Such 


behavior is expected in models where emission is pro¬ 
duced by curvature radiation (CR) occurring at the 
radiation-reaction limit in the outer magnetosphere 
Abdo et al.||2010b Harding et al.|2008 . This has be¬ 


come the most favored general description of gamma- 
ray emission from pulsars in the Fermi-hAT era. In 
these models one expects that the SED will fall off ex¬ 
ponentially above the break energy. There is, however, 
compelling evidence suggesting that CR occurring in 
the outer magnetosphere is not a complete description 
of pulsar emission at, and above, the GeV SED break: 

1. The discovery of power-law-type emission from 
the Crab pulsar at energies exceeding 100 GeV^ 
cannot be easily explained by curvature emis¬ 
sion from the outer magnetosphere [Aliu et al. 


2011 Lyutikov et al. 2012 unless the radius of 


curvature of the magnetic field line is larger than 
the radius of the light cylinder Bednarek|2012 


Some recent models attribute the pulsed very- 


^At this symposium the MAGIC collaboration presented ev¬ 
idence indicating that the power-law spectrum of the Crab pul¬ 
sar may extend to TeV energies. See http://fermi. gsfc .nasa. 
gov/science/mtgs/syiiiposia/2014/abstracts/185 


high-energy (VHE; E >100 GeV) emission from 
the Crab pulsar to inverse-Compton (IC) scat¬ 
tering originating in the outer magnetosphere 
[Du et al. 2012 Lyutikov et al. 2012[ Lyutikov 


2012 or to IC scattering from beyond the light 
cylinder Aharonian et al.||2012 Petri|[MT2 . 


2. The radiation-reaction limit of CR occurs when 
the acceleration gains achieved by an electron 
are equaled by radiation losses. The photon en¬ 
ergy at which this occurs in the outer magne¬ 
tosphere can be expressed in terms of the pul¬ 
sar period, the surface magnetic field strength, 
the radius of curvature of the accelerated par¬ 
ticle and an efficiency factor. [Lyutikov et al. 


2012 has shown that the break-energy values 


for several pulsars reported in the first Fermi- 
hAT pulsar catalog are so high that they require 
the efficiency factor and radius of curvature to 
approach or even reach their maximal allowable 
values^. 

3. Recent studies of the Geminga pulsar with 
Fermi-hAT and VERITAS (see Figure 1) show 
that the SED above the GeV break is compatible 
with a steep power law Aliu et al. 2015[ Lyu- 
tikov|2012 , but no emission has been seen above 
100 GeV. Similar conclusions can be drawn from 
an analysis of the Vela pulsar with Fermi-hAT 

who show that 


data from Leung et al. 2014 


multi-zone or time-dependent emission models 
are needed to fit the slower-than-exponential fall 
of the SED above 10 GeV. 

The question of whether the Grab pulsar is unique, 
or whether non-exponentially-suppressed gamma-ray 
spectra are common in gamma-ray pulsars is of great 
importance. Beyond the modeling of pulsar emission, 


^In more realistic models the acceleration efficiency is ex¬ 
pected to be a few percent to a few tens of percent |Lyutikov| 
|et al.|2012|. 
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Power-law x Exp. Cut-off fit [Likelihood] 
Power-law fit [Likelihood with E > 10 GeV] 
Power-law fit [Extrapolated above 100 GeV] 

Older Measurments of Geminga 
Singh et al, 2009 [PACT] 

Neshpor et al, 2001 [Crimea] 

Aharonian et al, 1999 [HEGRA] 

* Akerlof et al, 1993 [Whipple] 

Vishwanath et al, 1993 [Ootacamund] 

Crab Nebula and Pulsar 
— - - Crab Nebula [Albert et al, 2008] 

1% Crab Nebula 
Crab Pulsar [Aliu et al, 2011] 


(a) Phase-averaged 


(b) PI 


(c) P2 


[ht] 


Figure 1: SEDs data points and flux upper limits for the Geminga pulsar. Measurements of the Crab Nebula and 
pulsar are plotted for comparison. It is clear, even in the phase-resolved analysis, that the SED falls slower than an 


exponential and appears more consistent with a simple power-law. Figure taken from Aliu et al. 2015 


questions concerning the emission spectra of pulsars 
have significant implications for galactic dark matter 
searches, where unassociated gamma-ray excesses can 
be interpreted as the remnants of dark matter annihi¬ 
lation (e.g.,|Ab^ajian&^apIinghat]20^). Since pul¬ 
sars are likely the main background for these searches, 
categorizing the shape of pulsar spectra is a critical 
step towards validating any indirect dark matter sig¬ 
nal in the gamma-ray domain. To search for non- 
exponentially-suppressed emission above 50 GeV, we 
have performed a stacked analysis of gamma-ray pul¬ 
sars which is sensitive to emission which cannot be 
resolved in the Fermi-hAT analysis of individual ob¬ 
jects, but can be detected if aggregated from an en¬ 
semble. A stacked analysis which yields evidence of 
cumulative emission above 50 GeV would prove that 
some population of gamma-ray pulsars clearly exhibits 
non-exponentially-suppressed emission. This would 
indicate that inverse-Compton or wind-zone emission 
is common in gamma-ray pulsars and that pulsars con¬ 
tribute to the sub-TeV diffuse emission of the galaxy. 


2. Analysis 

2.1. The Aperture Photometry Method 

A maximum likelihood fitting procedure is typically 
employed when performing spectral analysis of Fermi- 
LAT data. The Fermz-LAT data can also be analyzed 
with an Aperture Photometry (AP) method where the 
raw event counts from a region of interest (ROI) are 
combined with a measure of the instrument exposure 
(cm^ s) to the region to determine the flux. This AP 
method is less sensitive and less accurate than the 
likelihood fitting procedure but it “provides a model 
independent measure of the flux” and it “is less com¬ 
putationally demanding”^. We demonstrate here that 
the AP method can be used to produce accurate SEDs 


^http://ferml.gsfc.nasa.gov/ssc/data/analysis/ 
scitools/aperture_photometry.html 
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Figure 2: Aperture photometry analysis steps for the Crab pulsar. Panel (a) plots the phase distribution (light curve) 
of the Crab pulsar from 5.2 years of Fermi-LAT observations. The Off phase range, [0.71 — 0.99], is defined in the 2nd 
Fermi-hNT catalog of gamma-ray pulsars (2PC). Panel (b) plots the distribution of photon energies for events which 
fell in the On and Off phase ranges. The Off events have been scaled by a which is the ratio of the On phase gate(s) 
size to the Off gate(s) size . Panel (c) shows the energy distribution of the excess events and panel (d) shows the 
significance of the excess in each energy bin. Panel (e) shows the Fermi-LAT exposure for the ROI used in each energy 
bin determined from gtexposure. In panel (f) the Crab pulsar AP SED is plotted alongside the Crab pulsar SED 
determined from a likelihood fit done in the 2PC. A broken power-law fit to Fermi-LAT and VERITAS data from I Aliu 


the AP flux to the 2PC flux in each bin, showing the level of agreement between the AP method and the likelihood 
method. This figure is taken from |McCami| [2014| . 


et al. 2011 is plotted, as well as the VERITAS >100 GeV bow-tie. Below the SED plotted in panel (f) is the ratio of 
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from multi-year pulsar data sets since an accurate de¬ 
termination of the background rate can be measured 
from the “Oj(f-pulse” phase range. The analysis pre¬ 
sented here proceeds as follows: 


1. Over the 100 MeV to 1 TeV energy range, 
logarithmically-spaced energy binning with 4 
bins per decade is chosen. 


2. An ROI is chosen around each pulsar with an 
energy-dependent radius. The radius chosen 
is three times the 68% point-spread-function 
(PSF) containment radius determined from a 
“front-conversion” Vela analysis by |Ackermann| 


et al.| [2013 . In order to maintain sufficient 


statistics at high energies, the radius of the ROI 
was fixed to 0.45° above 10 GeV. 


3. The Fermi-hAT analysis tools gtselect, 
gtmktime, gtbin and gtexposure are then run 
over each pulsar ROI for all observations per¬ 
formed within the period of validity of the pulsar 
timing solution. 


4. The photon event list is barycentered and phase- 


folded using the Tempo2 package Hobbs et al. 


2006 with the Fermi Tempo2 plugin and the cor¬ 


responding timing solution. 


5. Within each energy bin, a cut on phase is ap¬ 
plied and events which fall within the Off phase 
region and those which fall outside this region 
- the On phase region - are selected. The ratio 
of the size of the On phase range to the size of 
the Off phase range, defined as a, is then used 
to scale the number of event counts in the Off 
phase region (Nos) to the number in the On re¬ 
gion (Non). 


115 pulsars listed in the Second Aermi-LAT Catalog 
of Gamma-ray Pulsars Abdo et al.||2013 , which shall 
be referred to as 2PC throughout®. The 115 pulsar 
sample is composed of 39 millisecond pulsars and 76 
“young” non-recycled pulsars with an average data set 
spanning 4.2 yr®. The six AP analysis steps listed in 
Section [2T] were followed for each pulsar and using the 
resulting values of Non, Noff and T, it is quite simple 
to determine the total excess, 


Exjot = 

i=i 


( 1 ) 


the total exposure. 


1=1 


and thus, the average flux. 


( 2 ) 


Ex* 

Flux^^ = 

av /T-j 
• t.nt 


( 3 ) 


for N pulsars in a given energy bin, i. The significance 
of the total excess is determined by the generalized 
version of Equation 17 from Li & Ma 1983 (see Aha- 


ronian et al. 20041. In cases where the significance 
is less than 2a, the method of Helene 1983j is used 
to derive the 95% confidence-level upper limit on the 
total excess, which is in turn used to compute a flux 
upper limit. 


3. Results 


6. The number of excess pulsed events is then de¬ 
fined as Nex=Non—aNofi and the flux is Nex di¬ 
vided by the exposure (T) calculated in step 3 
using gtexposure. The significance of the ex- 
cess is calculated using Equation 17 fromjLi ^ 


Ma 1983 


Following this procedure one can derive the energy 
distributions for the On and Off phase regions, and 
the instrument exposure, for any pulsar. These dis¬ 
tributions and the derived AP SED are shown for the 
Crab pulsar in Figure 2. 


The stacking analysis results for the young pul¬ 
sar and millisecond pulsar ensembles are shown in 
Figure 3. No significant excesses are seen in these 
analyses at energies above 50 GeV. Upper limits on 
the average flux, determined at the 95% confidence- 
level, are listed in Table |T| for three energy bins above 
50 GeV. Limits are also presented in units of the 
Crab pulsar where the broken power-law fit to the 
Fermi-LAT and VERITAS data presented in [Aliuj 
et al. 2011 defines a Crab pulsar unit. In addi¬ 


tion to these analyses, we stacked sub-samples of the 


2.2. Stacking The Pulsar Data Sets 

Fermi-hAT has detected over 150 new gamma-ray 
pulsars^ and the stacking performed in this work uses 


https: //confluence. slac. stcuif ord. edu/display/ 
GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars 


total of 117 pulsars are listed in the 2PC, however, the 
Crab pulsar was excluded from this analysis since we are inves¬ 
tigating whether high-energy Crab-pulsar-like emission is seen 
in other pulsars. Further, PSRJ2215-I-5135 was also excluded 
from the study since no Off phase region was listed for this 
source in the 2PC. 

®The amount of data analyzed here depends entirely on the 
availability and validity of pulsar spin-down timing solutions 
used for phase-folding. 
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Figure 3: Panel (a) shows the average flux (square markers) from 76 young pulsars determined by dividing the total 
excess by the total exposure (see Equations 1—3). The dashed-line histogram shows one over the total exposure, 
indicating the flux which would correspond to a single excess photon. This is the minimum possible flux which could be 
measured given the total exposure. The gray cross shows the most constraining limit on emission from a single pulsar 
in the 56.2—100 GeV range presented in the 2PC. The 2PC presented no limits at higher energies. The broken 
power-law fit to the Crab pulsar data from Aliu et al. [2011 is plotted for scale. Panel (b) plo ts the same qua ntities for 


the stacked analysis of 39 millisecond pulsars. This figure is adapted from figures presented in McCann 2014 


data where each sub-sample was composed of the 10 
pulsars with the largest value of a given parameter. 
Sub-sample selections based on gamma-ray luminos¬ 
ity, spin-down power, spin-down power over distance 
squared, gamma-ray photon flux and non-thermal X- 
ray energy flux were investigated^. No significant ex¬ 
cesses were observed above 50 GeV in any of these 
sub-sample stacking analyses. 

The shape of the average young pulsar and average 
millisecond pulsar SEDs were categorized by fitting a 
power law times a super-exponential cutoff function 


■>dF 

= A 
dE 


E 

1 GeV 


r 

-( E 
g V -^cut 




(4) 


to the SED data. These fits are presented in Fig¬ 
ure 1^ Fixing 6 = 1 reduces Equation 4 to a power 
law times an exponential cutoff function and, as ex¬ 
pected, this functional form does not reproduce the 
sub-exponential fall of the SED above the break. How¬ 
ever it can be used to measure the average flux- 
weighted value of the spectral index (T) and cutoff 
(Ecut) parameters Abdo et al.||2013| . It is clear from 
Figure]^ that the average SEDs have qualitatively the 
same shape, with the average flux from the 39 mil¬ 
lisecond pulsars about an order of magnitude lower 
than the average flux from the 76 young pulsars. The 
spectral parameters derived from fitting the two en- 


^The Crab pulsar was excluded from all of these sub-sample 
stacking analyses. The parameter values listed in the 2PC cat¬ 
alog were used in all cases. 


sembles are remarkably similar and are presented in 
the caption of Figure]^ 


4. Discussion and Conciusion 

Following a stacked analysis of 115 gamma-ray pul¬ 
sars, with an average exposure of ^4.2 yr per pul¬ 
sar, we find no evidence of cumulative emission above 
50 GeV. Stacked searches exclusive to the young pul¬ 
sars, the millisecond pulsars, and several other promis¬ 
ing sub-samples also return no significant excesses 
above 50 GeV. Any average emission present in the 
entire pulsar sample is limited to be below ^7% of 
the Crab pulsar in the 56-100 GeV band. The av¬ 
erage flux limits presented in Table |l] are roughly 3 
times lower than the best flux limits achieved in ded¬ 
icated individual pulsar analyses done in the 2PC in 
the 56-100 GeV band. 

One should note that a limit on the average flux 
from 115 pulsars at 7% of the Crab pulsar level is 
consistent with, for example, a scenario in which all 
115 pulsars emit at 7% of the Crab pulsar level. It 
is also consistent with a scenario in which 8 pulsars 
emit at 100% the level of the Crab pulsar and the 
remaining 107 pulsars have zero emission. Therefore 
this analysis does not exclude the possibility of finding 
several pulsars which are as bright as the Crab pulsar 
above 50 GeV, or several dozen which are ten times 
dimmer. It does, however, constrain the average flux 
from the ensemble, and therefore for every individ¬ 
ual pulsar detected above this flux limit, the average 
emission from the remaining pulsars is constrained to 
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All 

Young Pulsars 

Millisecond Pulsars 

Energy Range 

Flux Limit 

Flux Limit 

Flux Limit 

Flux Limit 

Flux Limit 

Flux Limit 

[GeV] 

[xlO"^^ 

[Crab pulsar 

[xlO'^^ 

[Crab pulsar 

[xlO"^^ 

[Crab pulsar 



units] 

cm“^s“^] 

units] 

cm“^s“^] 

units] 

56.2 — 100 

1.57 

0.07 

2.03 

0.09 

1.44 

0.07 

100 — 177 

1.52 

0.31 

1.88 

0.38 

1.14 

0.23 

177 — 316 

1.34 

1.21 

1.96 

1.76 

0.50 

0.45 


Table I Limits at the 95% confidence level on the average flux from stacked ensembles of gamma-ray pulsars. The limit 


values presented in Crab pulsar units assume the broken power-law fit to the Crab pulsar data from Aliu et al. 
a Crab pulsar flux unit. 


2011 



Energy [MeV] 


Figure 4: The average SEDs derived from the stacking of 
the 76 young pulsars and 39 millisecond pulsars. The 
SEDs are each fit with a power law times a 
super-exponential cutoff keeping b both fixed to unity 
and allowing it to float. For the pure exponential cutoff 
case {b = 1) the best fit F value is 0.54T0.05 for the 
millisecond pulsars and 0.41T0.01 for the young pulsars 
while the best fit Ecut values are 3.60±0.21 GeV and 
3.54T0.04 GeV, respectively. Allowing b to float we find 
that sub-exponential forms (6 < 1) are preferred, with 
the best-fit b value of 0.59±0.02 for the young pulsars 
and 0.7±0.15 for the millisecond pulsars. The broken 
power-law fit to the Crab pulsar data from | Aliu et al.| 
s is plotted for scale. Note that only statistical 
uncertainties on the SED data points were used during 
the fitting and thus the uncertainty on the best-fit 
parameter values are likely underestimated. This figure is 
taken from McCann [ 20 T 4 ] . 


be further below the limit. 

In the 100 MeV to ~50 GeV energy range we find 
that the average SEDs returned from the young pul¬ 
sar and millisecond pulsar stacking analyses are very 
similar in shape and are generally compatible with a 
power law times a sub-exponential cutoff. |Abdo et al.| 


2010 and Celik & Johnson 2011 have shown that 


a sub-exponential cutoff function approximates a su¬ 
perposition of exponential cutoffs, thus the appear¬ 
ance of a sub-exponential cutoff in the ensemble SED 


is to be expected within a curvature radiation model. 
We note, however, that the highest energy spectral 
point is higher than the best fit sub-exponential cutoff 
function at the ~2.4cr level in both the young pulsar 
and millisecond pulsar cases. This cannot be taken 
as strong evidence for a non-exponentially-suppressed 
pulsar emission component aggregating in the stacked 
analysis, however, the available data cannot rule it out 
beyond the level of the limits shown in Figure 3 and 
Table m 

Beyond this work, improvements can be made using 
the forthcoming Fermi-LAT pass-8 data release which 
will improve the Fermi-hAT acceptance by ^25% at 
100 GeV Atwood et al.| 2013 . Improvements to this 
stacking analysis can also be made by employing a 


likelihood framework to stack the sources (see Acker- 


mann et al. 2011 for example), rather than the sim¬ 


ple On minus Oj(f procedure described here. The flux 
sensitivity of any stacking analysis will, however, ul¬ 
timately be bounded by the exposure of the Fermi- 
LAT. A future stacking analysis which doubles both 
the number of pulsars and the duration of observa¬ 
tion used will increase the exposure by a factor of 4, 
indicating that future stacking analyses which do not 
yield detections may improve on the limits presented 
here by perhaps one or two orders of magnitude. 

A more detailed account of the stacking analysis 
methods and results of this study can be found in 


McGann 2014 
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